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Abstract 

Due to their inherent dissipation and stability, the MacCormack scheme and its variants 
have been widely used in the computation of unsteady flow and acoustic problems. 
However, these schemes require many points per wavelength in order to propagate waves 
with a reasonable amount of accuracy. In this work, the linear wave propagation 
characteristics of MacCormack-type schemes are investigated, and methods for greatly 
improving their performance are described and demonstrated. 


Introduction 

In the field of computational aeroacou sties, numerical schemes are expected to propagate 
waves accurately for long distances over long periods of time. In order to accomplish this 
goal, a certain number of spatial points are required per wavelength to model each wave, 
and a certain time step is required in order to model the wave’s movement in time. It is 
desirable from a computational standpoint to reduce the number of points required per 

wavelength and increase the size of the allowable time step. 



One popular and well-tested method uses a modification of the MacCormack scheme [1], 
which is second order accurate in time and fourth order accurate in space. This extension 
of the MacCormack scheme is known as the 2-4 scheme, and was described by Gottlieb 
and Turkel [2]. This scheme has been used successfully on a wide range of fluid and 
aeroacoustics problems [3-15]. Sankar, Reddy, and Hariharan have evaluated this scheme 
for aeroacoustics applications [16]. It has been extended to sixth- order spatial accuracy 
by Bayliss, et. al. (2-6 scheme) [17]. 

In this paper, these MacCormack-type schemes will be investigated in detail, and 
extensions to high-order accuracy will be developed. A test problem is used to quantify the 
performance of the various schemes. These schemes have been validated on the real-world 
problem of noise radiated by a supersonic jet using the linearized Euler equations. 


Test Equation and Numerical Formulation 

For comparative purposes, a simple hyperbolic equation is solved: 

dU_ = _ dU_ 
dt dx 


(1) 


The harmonic solution to this equation is: 

( 2 ) 

where k is the wavenumber and co is the frequency. Here, 

w = ck. (3) 

In order to numerically solve Eq. (1), the equation must be discretized in time and space. 
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The discretized spatial derivatives must model accurately the derivative of the waveform: 

(4) 


— = iAe i( “ ) = — YajUj 
dx A - ^ J J 


AX ; 


j=-N 


where the coefficients aj are the coefficients used in taking the spatial derivative. 


Since: 


jj _ e i{k{x+jbx)-at) 


__ gikj&x 


(5) 


numerical wavenumber can be defined, following Tam and Webbs work [18]: 


• M 

l V' ifc/Ax 

Ax ^ * 


= — y\ a i e 

Av ^ J 


( 6 ) 


In this way, Eq. (4) becomes: 


— = ik V (fa_rt) 


(7) 


As the numerical scheme marches in time, the time integration must also model accurately 
the evolution of the waveform. From before: 


dU _ <?£/ 

dt dx 

= -icke i{kx -°* ) 

= - iox 


( 8 ) 


Integrating Eq. (8) gives: 


!/(*,« + AO = e^V (fat " -) 

=u(x,ty i[0 * t) 


(9) 


One popular method for numerical time integration is the Runge-Kutta scheme. A generic 
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six-stage Runge-Kutta method has this form: 

U ( 1 ) =U(x,t) 
dU w 


U m = U(x,t) - a 2 cAt 
U m = U(x,t) - a 3 cAt 
U (4) = U(x,t) - a 4 cAt 
U (S) = U{x,t)-a s cAt 
U i6) = U(x,t)-a 6 cAt 


dx 

dlf^_ 

dx 

dV r(3) 
dx 

dV r(4) 
dx 

dl r (5) 

dx 


U(x,t + At) = U(x,t)~ cAt 


dU^_ l) . „ dU w „ dU w 
dx 


n w w o ^ n 


dx 

a dU w 0 dU m a dU (6) 

+ A— + A— +A— j 


( 10 ) 


Using the coefficients from Table I for the sixth-order accurate Runge-Kutta scheme and 
using Eq. (7) to define the numerical wavenumber, Eq. (10) becomes: 


U(x,t + At) = U(x,t) 


, (-ick* Atf ( -ick* At ) 3 

1 - ick At + i i —- 1 - 

2 6 

(-ick* At) (— ick* At) (-ick* At) 

+ - — + - — + - — 

24 120 720 


(ID 


= U(x,t)e~ i6} ' M 


Comparing Eq. (1 1) to Eq. (9), we can define the numerical frequency co* as: 


(O* 


-i . U(x,t + At) 
In . r 
At U(x,t) 


( 12 ) 
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In the MacCormack scheme, integration in time is accomplished by applying an operator to 
the solution at the last rime level. The time integration method used in the existing schemes 
is a second-order accurate Runge-Kutta method: 


U (1) = U(x,t) 

U (2) =U(x,t)-cAt- 


At 

'dU (1) 

2 

dx 


In the MacCormack scheme, one-sided differences are used in order to add dissipation to 
stabilize the scheme. The one-sided differences are defined in this way: 

dU F & TT 

■si 

( 14 ) 

dU B ^ TT 

*■, -k a ' u - j 

such that the underlying central difference is recovered when the forward and backward 
differences are added together: 


f f ■ iw. 

UX i j=-M 

Aj=±(aj+ a _j) 


For example, in the original MacCormack scheme, forward and backward differences are 
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defined as first-order accurate differences: 


f- =jz(- u ‘ + u ^ 

dx i Ax 

dU . <? 2 £7 2 x 

= ar +A *l^ +,?(Al) 

dx i Ax 

dU A d 2 U Q / A 2 V 

‘-ar^sr*** > 


It is seen that adding (16) and (17) together recovers a second-order accurate central 
difference. In all of the existing MacCormack-type schemes, the one-sided differences are 
first-order accurate, with a dissipative leading enror. 

The forward and backward differences are alternately used in the time-stepping routine. 


For example, Eq. (13) becomes: 

U (1) = U{x,t) 

U (2) =U{x,t)-cAt- 


U(x,t + At) = U{x,t)-c 


At f dU (1) F dU m 8 


2 dx 
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Using Eq. (6), the numerical wavenumbers can be found for the one-sided differences as 
well as the underlying central difference: 

k' F =^(e^-l) = (k' + i6) 

k l = j^{- e ~ ik6X + 1 ) = ( k '- i5 ) (19) 

kc = —(e ik ^-e- ikAx )^k' 

c Ax v ' 

Notice that the real part of the numerical wavenumber is identical to that of the underlying 
central difference for all differences, and that the imaginary pans are due to the form of the 
one-sided differences and are equal and opposite. 

Putting these definitions into the time integration method gives: 

U (l] =U(x,t) 

U i2) = U(x,t) - icAtk F U (1) (20) 

= (l-icA tk F ^)U(x,t) 

Then: 

U(x,t + At) = U(x,t) - icf[k' F U w + k' B U w ] 

= U{x,t)-ic^-[k;U{x,t) + k' B ( l-icAtk' F )U(x,t)\ 

= i -ic^(k;+ki)-^^(kik;)u(x,t) (2D 

= 1- tcA^‘ + ( ~^ )2 (Jfe*) 2 + 5 2 ) u{x,t) 
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Using Eq. (13), it is clear that the numerical frequency for the MacCormack scheme is: 


(o' 



1 - icAtk' 



( 22 ) 


From this development, three contributions to the error are present in the MacCormack 
scheme: the error in dispersion from the spatial central difference, the error in dissipation 
from the one-sided spatial differences, and the error in both dispersion and dissipation from 
the time integration method itself. 


Higher Order Runge-Kutta Schemes for MacCormack-tvpe Methods 

In order to investigate the effect of the time integration scheme on the dispersion and 
dissipation error, several high-order Runge-Kutta methods were constructed for 
MacCormack-type schemes. Eq. (10) gives the generic form of the Runge-Kutta schemes. 
For example, the two-stage, second-order accurate method given in Eq. (13) would have 
the coefficients given for the RK2 scheme. 

It should be noted that, in order to achieve the desired spatial accuracy using a 
MacCormack-type scheme, the sum of the odd P's (1,3,5) must add up to equal the sum of 
the even p's (2,4,6). In this way, the central difference will be obtained from the sum of 
the forward and backward differences. Similar constraints are used for the odd time 
derivatives to insure that the one-sided differences sum properly. 

Table 1 gives the coefficients for Runge-Kutta 2nd, 4th, and 6th order accurate methods, as 
well as the coefficients for Hu, et. al.'s optimized 4-6 Runge-Kutta scheme [19]. 
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Table I: Coefficients for Runge-Kutta Time Stepping Schemes. 

Following the development given for the 2-4 scheme, we find that a sixth order accurate 
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Runge Kutta method gives this result: 


l-ic&tk‘ + ^ lC ^ (Jfe*) 2 +5 2 j + 
w * = Zii n J k'((k'f + S 2 )+ -l^ ---{(k') 2 + S 2 ) 2 + (23) 

Again, it is seen that the errors arise from the dispersion error in the central difference and 
the dissipation error in the one-sided differences. 

From Eq. (23), it can be seen that the numerical frequency for a MacCormack-type scheme 
may be written generally as: 



where S is the number of stages in the Runge-Kutta scheme and the ci’s are the leading 
coefficients in the time integration, given in Table I. 


Improved Central Differences 


Various central differences were investigated during this work for use with MacCormack- 
type schemes. Fourth, sixth, and eighth order central differences were investigated, as 
well as the fourth-order accurate Dispersion Relation Preserving (DRP) scheme of Tam and 
Webb. 
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The error in the dispersion relations was defined as: 

Error - k* Ax - kAx (25) 

The dispersion errors are plotted as a function of wavenumber kAx in Figures 1 and 2. It 
can be seen that higher orders of accuracy result in lower errors, as would be expected. 
However, lower order differences can be optimized as shown by Tam and Webb to provide 
better performance for a given level of acceptable error, illustrated in Figure 2. 


Improved One-Sided Diff erences 

1 Higher-Order Accurate One-Sided Differences 

Using the basic methodology of the MacCormack method, improved accuracy can be 
attained by altering the one-sided differences to have higher accuracy while still adding up 
to a more accurate central difference. For a fourth-order accurate central difference, two 
schemes can be constructed: 

The first scheme is the normal 4/2 scheme, which is fourth order accurate in space and has 
second-order dissipation in each step. The next scheme is a 4/4 scheme, which is fourth 
order accurate in space and has fourth-order dissipation in each step. This higher-order 
accuracy is achieved by adding one point to each of the one-sided differences, as shown in 

Table II. 

Notice how both schemes add up to the identical 4th order accurate central difference. The 
difference in the two schemes is in the leading error terms of the one-sided differences, 
which affects the inherent dissipation of the scheme. The effect of this change is to lower 
the inherent dissipation at a given wavenumber, as shown in Figure 3. This gives the 
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O-i 


a, 

a 2 

a 3 

a 4 




-1 

1 


0 

0 

2/2 [1] 

0 

0 

Ax 

Ax 

0 




-7 

8 

-1 

0 

0 

4/2 [2] 

0 

0 

6 Ax 

6 Ax 

6 Ax 

6/2 [17] 

0 

0 

-37 

45 

-9 

1 

0 

30 Ax 

30 Ax 

30 Ax 

30 Ax 


0 

0 

-533 

672 

-168 

32 

-3 

8/2 



420Ax 

420Ax 

420Ax 

420Ax 

420Ax 



-2 

-3 

6 

-1 

0 

0 

4/4 

0 

6 Ax 

6 Ax 

6 Ax 

6 Ax 



-9 

-19 

36 

-9 

1 

0 

6/4 

0 

30 Ax 

30 Ax 

30 Ax 

30 Ax 

30 Ax 



-.3766 

-.4968 

11651 

-.3334 

.04168 

0 

DRP/4 

0 

Ax 

Ax 

Ax 

Ax 

Ax 



-.30874 

-.6326 

L2330 

-.3334 

.04168 

0 

DRP/opt 

0 

Ax 

Ax 

Ax 

Ax 

Ax 



-120 

-293 

552 

-168 

32 

-3 

8/4 

0 

420Ax 

420Ax 

420Ax 

420Ax 

420 Ax 

420Ax 


3 

-30 

-20 

60 

-15 

2 

0 

6/6 

60 Ax 

60 Ax 

60 Ax 

60 Ax 

60 Ax 

60 Ax 


18 

-192 

-185 

480 

-150 

32 

-3 

8/6 

420Ax 

420Ax 

420Ax 

420Ax 

420Ax 

420Ax 

420Ax | 


Table II: Coefficients for MacCormack-Type Schemes 


scheme a wider range of wavenumbers that it can accurately resolve, and thus requires less 
points per wavelength. 
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Using this basic methodology, families of 6th and 8th order accurate schemes can be 
constructed. In these higher order schemes, it is possible to add another point to the one- 
sided differences to attain 6th order dissipation. These results are given in Table II. 

2 . Optimization of the Split Operators 

The Tam and Webb optimization technique may be used on the one-sided differences to 
improve their performance at wavenumbers of interest. The procedure is as follows: 

The central difference that the one-sided differences must add up to is fixed: 


dU 

c 

A_ 3 L r i _ 3 + A_ 2 U t _ 2 + A.!?/,.! + A 0 U i 

dx 

i 

+ Ajt/j +1 + AqUi + 2 "h A3L/ "j +3 


This, in turn, determines some of the coefficients of the one-sided difference: 

= 2 A 3 U i+ 3 + 2 A 2 U i+ 2 +a 1 U i+l + a 0 U i + a_ l U i _ l ( 27 ) 

dx : 


There are two conditions that are known for the three unknowns: 

2 A3 + 2 Aj + a x + a 0 + a_ t = 0 

a 1 -a_i = 2 Aj 


( 28 ) 


For the third condition, we can set a range to minimize the dissipation of the one-sided 
differences. For this work, the dissipation was minimized at 8 points per wavelength and 
higher. To do this, the following equation is minimized (following Tam and Webb's 
procedure): 


E=J 


4 


3 1 

i 

Re 


TC 




die 


4 


( 29 ) 
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The condition that E is a minimum is: 


dE 
da , 


= 0 


(30) 


This gives the condition required, and the coefficients are solved from the resulting linear 
system of equations. Figure 3 illustrates the effect of optimization on the inherent 
dissipation of the scheme. Notice that there is a price for this optimization: the lower 
wavenumbers (more points per wavelength) have slightly increased dissipation as 
compared to the original scheme. 

The DRP schemes shown are given in Table II. 

3. Dissipative Errors of the One-Sided Differences 

Figure 3 shows the dissipative errors in the one-sided differences for various 
MacCormack-type schemes. It can be seen that, in general, the dissipative performance of 
schemes with 2nd order dissipation (4/2 and 6/2, for example) are similar. It can also be 
seen that the closer the split stencil is to the central difference stencil, the less dissipation the 
scheme will have (6/6, for example). It can be seen that, in general, higher-order accuracy 
reduces dissipation errors at a given wavenumber. The positive effect of the DRP-style 
optimization on the dissipation errors can also be seen in this figure. 


Performance of the Improved MacCormack-tvpe schemes 

In this section, the total performance of each scheme will be quantified. In order to fairly 
compare the schemes, the errors in dissipation and dispersion per wavelength of travel will 
be shown as functions of wavenumber and time step. In this way, the effects of both the 
time and space discretization are shown. 
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In these results, the following definitions are used: 


PPW = 


2k 

kAx 


CFL = 


cAt 


NST = 


Ax 

PPW 

CFL 


( 31 ) 


Here, PPW is the number of spatial points per wavelength, CFL is the Courant number, 
and NST is the number of time steps per cycle. 


From Eq. (1 1), we see that: 

U{x,t + NST At) = U{x,ty i{o>NS ™ ] 

= U(x4e-^) NST 


(32) 


Since the problem is periodic, 

U(x,t + NSTAt) = U(x,t) 

(33) 

(e _ ‘ ( ^)) = 1 + Oi 

and we can define the numerical solution to be: 

= A +iD (34) 

Here, A is the amplitude per wavelength traveled and D is the dispersion error per 
wavelength traveled. 

Figure 4 shows the amplitude and dispersion per wavelength traveled for the previously 
published MacCormack-type schemes. It is seen that there is a large amount of dissipation 
inherent in these schemes, acting over a large part of the wavenumber range. Also, the 
dispersion plots emphasize the low order of accuracy in time. 
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Figure 5 illustrates the effect that varying the time integration method has on the 
MacCormack-type schemes. In Figure 5, the 6/4 one-sided differences are used, and 
various time integration methods are employed. Notice how the variation of the errors with 
increasing time step is greatly reduced with the 46 and 6 stage time integration methods. 
Also, the optimized 46 time integration method is very nearly as accurate as the 6-stage, 6th 
order method while requiring less work. 

Figure 6 compares the 2-4/2 method with the best of the two newer methods: the 46-6/4 
method and the 46-DRP/optimized method. Notice that using the DRP method can extend 
the usable wavenumber range for a given amount of error significantly without requiring 
any additional calculations. Also, the stability range of the newer schemes are much greater 
than that of the older MacCormack-type schemes, while retaining very good accuracy. 

Figures 7-9 shows the final results of this work: the actual wave propagation. 

In Figure 7, a wave of 8 points per wavelength (kAx = .785) is propagated for 400 time 
steps at a CFL of 0.5 (25 wavelengths of travel). Results are shown for the original 2-4/2 
scheme, the 46-6/4 scheme, and the optimized 46-DRP/opt scheme. The low CFL number 
is chosen due to the low stability limits of the 2-4/2 scheme. 

It can be seen that the 2-4/2 scheme is the most dissipative, followed by the 46-6/4 scheme 
and the 46-DRP/opt scheme. Referring to Figure 3, the magnitudes of the dissipation for 
these schemes fall in the same sequence.. 

The dispersive error is illustrated in Figure 1. The fourth-order central difference has the 
highest dispersive error, with the numerical wave trailing the actual wave. The sixth-order 
central difference and the DRP central difference have about the same amount of 
dispersion, with the sixth order wave trailing the actual wave, and the DRP wave leading. 
Again, Figure 1 illustrates this effect. 
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Moving on to higher CFL, Figure 8 shows the same problem with a CFL of 1.25 (8 points 
per wavelength, 25 wavelengths of travel, 160 time steps). The 2-4/2 scheme was unstable 
at this larger time step, so only the results from the other two schemes are shown. It can be 
seen that the results are almost identical to the previous calculation, which illustrates the 
large time steps possible with the optimized 46 time stepping scheme. 

Figure 9 shows the results for 6 points per wavelength (kAx = 1.05) after 33.33 
wavelengths of travel (CFL = 1.25, 160 time steps). The improvement in both dispersion 
and dissipation of the DRP scheme can be seen compared to the 46-6/4 scheme, illustrating 
the benefits of optimization for this severe problem. 


Conclusions 

In this work, the dispersive and dissipative characteristics of the existing MacCormack-type 
scheme were investigated, and several ways were found to improve the accuracy of this 
type of scheme. The accuracy and stability of these schemes were greatly enhanced 
without much additional calculations being needed. 

The MacCormack-type schemes are of great interest due to their ease of programming and 
use, and inherent numerical dissipation. This work shows that this type of scheme can be 
optimized to perform very well. 
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k* ax - kAx 



Figure 1 

Effect of Central Difference on the Onset of Dispersion Error. 




|kAx - k*Ax| 
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Figure 2 

Effect of Central Difference on the Magnitude of the 

Dissipation Error. 
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Figure 3 

Effect of One-Sided Differences on the Magnitude of the 

Dissipation Error. 
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Figure 4 

Performance of Previous MacCormack-Type Schemes. 


23 



kAx 



Figure 5 

Effect of Time Integration Method on Accuracy. 
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Figure 6 

Comparison of Accuracy of MacCormack-Type Schemes. 
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Figure 7 

Comparison of results obtained by MacCormack-type methods 
(8 ppw, CFL = 0.5, 25 wavelengths of travel) 
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Figure 8 

Comparison of results obtained by MacCormack-type methods 
(8 ppw, CFL = 1.25, 25 wavelengths of travel) 
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Figure 9 

Comparison of results obtained by MacCormack-type methods 
(6 ppw, CFL = 1.25, 33.3 wavelengths of travel) 
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